function s = smoothdata3(d, f, t)
% smoothdata3  - apply gaussian smoothing to 3D data
%
% FORMAT:       s = smoothdata3(d, f [, t])
%
% Input fields:
%
%       d           3D double data (for N-dim, each 3D volume is smoothed)
%       f           FWHM kernel sizes
%       t           smoothing kernel threshold (default: 0)
%
% Output fields:
%
%       s           smoothed data

% Version:  v0.7g
% Build:    9031300
% Date:     Mar-13 2009, 12:00 AM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
   (~isa(d, 'double') && ...
    ~isa(d, 'single')) || ...
    ndims(d) < 3 || ...
   ~isa(f, 'double') || ...
    numel(f) ~= 3 || ...
    any(isinf(f) | isnan(f) | f < 0 | f > 32)
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing argument.' ...
    );
end
if nargin < 3 || ...
   ~isa(t, 'double') || ...
    numel(t) ~= 1 || ...
    isinf(t) || ...
    isnan(t) || ...
    t < 0
    t = 0;
elseif t > 0.25
    t = 0.25;
end
f = f(:)';

% build 3 separate smoothing kernels for each dim
k1 = smoothkern(f(1), t);
k2 = smoothkern(f(2), t);
k3 = smoothkern(f(3), t);

% build 3D kernels
kc1 = zeros(max(max(numel(k1), numel(k2)), numel(k3)) * ones(1, 3));
kc2 = kc1;
kc3 = kc1;

% fill kernels
nk = size(kc1, 1);
mk = (nk + 1) / 2;
kc1(1+(nk-numel(k1))/2:(nk+numel(k1))/2, mk, mk) = k1;
kc2(mk, 1+(nk-numel(k2))/2:(nk+numel(k2))/2, mk) = k2;
kc3(mk, mk, 1+(nk-numel(k3))/2:(nk+numel(k3))/2) = k3;

% apply as three steps, but in one go
if ndims(d) == 3
    s = conv3d(conv3d(conv3d(d, kc1), kc2), kc3);
else
    s = d;
    sz = size(s);
    ns = prod(sz(4:end));
    for sc = 1:ns
        s(:, :, :, sc) = ...
            conv3d(conv3d(conv3d(s(:, :, :, sc), kc1), kc2), kc3);
    end
end
